Brian E. J. Rose, University at Albany
You really should be looking at The Climate Laboratory book by Brian Rose, where all the same content (and more!) is kept up to date.
Here you are likely to find broken links and broken code.
This document uses the interactive Jupyter notebook format. The notes can be accessed in several different ways:
github at https://github.com/brian-rose/ClimateModeling_coursewareAlso here is a legacy version from 2015.
Many of these notes make use of the climlab package, available at https://github.com/brian-rose/climlab
In [1]:
# Ensure compatibility with Python 2 and 3
from __future__ import print_function, division
Suppose that CO$_2$ actually behaved as a grey gas. In other words, no spectral dependence in absorptivity.
If we then double the CO2 concentration in the atmosphere, we double the number of absorbers. This should imply that we also double the absorption cross-section:
$$ \kappa^\prime = 2 ~ \kappa $$This would imply that we double the optical thickness of every layer:
$$ \Delta \tau^\prime = 2 \left( -\frac{\kappa}{g} \Delta p \right) = 2 ~ \Delta \tau$$And since (from Lecture 9) the absorptivity / emissivity of each layer is
$$ \epsilon = 1 - \exp\big( - \Delta \tau \big) $$the modified absorptivity is
$$ \epsilon^\prime = 1 - \exp\big( - 2\Delta \tau \big) = 1 - \left( \exp\big( - \Delta \tau \big)\right)^2 = 1 - (1-\epsilon)^2 $$or simply $$ \epsilon^\prime = 2 \epsilon - \epsilon^2 $$
(Note that $\epsilon^\prime = 2 \epsilon$ for very thin layers, for which $\epsilon$ is small).
In [2]:
# Applying the above formula
eps = 0.586
print( 'Doubling a grey gas absorber would \
change the absorptivity from {:.3} \
to {:.3}'.format(eps, 2*eps - eps**2))
If CO2 behaved like a grey gas, doubling it would cause a huge increase in the absorptivity of each layer!
Back in Lecture 7 we worked out that the radiative forcing in this model (with the observed lapse rate) is about +2.2 W m$^{-2}$ for an increase of 0.01 in $\epsilon$.
This means that our hypothetical doubling of "grey CO$_2$" should yield a radiative forcing of 53.5 W m$^{-2}$.
This is an absolutely enormous number. Assuming a net climate feedback of -1.3 W m$^{-2}$ K$^{-1}$ (consistent with the AR5 ensemble) would then give us an equilibrium climate sensitivity of 41 K.
This figure shows the solar radiation spectrum for direct light at both the top of the Earth's atmosphere and at sea level. The sun produces light with a distribution similar to what would be expected from a 5525 K (5250 °C) blackbody, which is approximately the sun's surface temperature. As light passes through the atmosphere, some is absorbed by gases with specific absorption bands. Additional light is redistributed by Raleigh scattering, which is responsible for the atmosphere's blue color. These curves are based on the American Society for Testing and Materials (ASTM) Terrestrial Reference Spectra, which are standards adopted by the photovoltaics industry to ensure consistent test conditions and are similar to the light that could be expected in North America. Regions for ultraviolet, visible and infrared light are indicated.
Source: http://commons.wikimedia.org/wiki/File:Solar_spectrum_en.svg
Careful: I'm pretty sure what is plotted here is not the total observed spectrum, but rather the part of the emissions from the surface that actual make it out to space.
As we now, the terrestrial beam from the surface is depleted by absorption by many greenhouse gases, but principally CO$_2$ and H$_2$O.
However there is a spectral band centered on 10 $\mu$m in which the greenhouse effect is very weak. This is the so-called window region in the spectrum.
Since absorption is so strong across most of the rest of the infrared spectrum, this window region is a key determinant of the overall greenhouse effect.
Another big shortcoming of the Grey Gas model is that it cannot represent the water vapor feedback.
We have seen above that H$_2$O is an important absorber in both longwave and shortwave spectra.
We also know that the water vapor load in the atmosphere increases as the climate warms. The primary reason is that the saturation vapor pressure increases strongly with temperature.
In [3]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import xarray as xr
from xarray.ufuncs import cos, deg2rad, log
# Disable interactive plotting (use explicit display calls to show figures)
plt.ioff()
In [4]:
# Open handles to the data files
datapath = "http://thredds.atmos.albany.edu:8080/thredds/dodsC/cesm/"
ctrl = xr.open_dataset(datapath + 'som_control/som_control.cam.h0.clim.nc', decode_times=False)
co2 = xr.open_dataset(datapath + 'som_2xCO2/som_2xCO2.cam.h0.clim.nc', decode_times=False)
In [5]:
# Plot cross-sections of the following anomalies under 2xCO2:
# - Temperature
# - Specific humidity
# - Relative humidity
fig, axes = plt.subplots(1,3, figsize=(16,6))
ax = axes[0]
CS = ax.contourf(ctrl.lat, ctrl.lev, (co2['T'] - ctrl['T']).mean(dim=('time','lon')),
levels=np.arange(-11,12,1), cmap=plt.cm.seismic)
ax.set_title('Temperature (K)')
fig.colorbar(CS, orientation='horizontal', ax=ax)
ax = axes[1]
CS = ax.contourf(ctrl.lat, ctrl.lev, (co2['Q'] - ctrl['Q']).mean(dim=('time','lon'))*1000,
levels=np.arange(-3,3.25,0.25), cmap=plt.cm.seismic)
ax.set_title('Specific humidity (g/kg)')
fig.colorbar(CS, orientation='horizontal', ax=ax)
ax = axes[2]
CS = ax.contourf(ctrl.lat, ctrl.lev, (co2['RELHUM'] - ctrl['RELHUM']).mean(dim=('time','lon')),
levels=np.arange(-11,12,1), cmap=plt.cm.seismic)
ax.set_title('Relative humidity (%)')
fig.colorbar(CS, orientation='horizontal', ax=ax)
for ax in axes:
ax.invert_yaxis()
ax.set_xticks([-90, -60, -30, 0, 30, 60, 90]);
ax.set_xlabel('Latitude')
ax.set_ylabel('Pressure')
fig.suptitle('Anomalies for 2xCO2 in CESM slab ocean simulations', fontsize=16);
In [6]:
fig
Out[6]:
In fact the specific humidity anomaly has roughly the same shape of the specific humidity field itself -- it is largest where the temperature is highest. This is a consequence of the Clausius-Clapeyron relation.
The relative humidity anomaly is
The smallness of the relative humidity change is a rather remarkable result.
This is not something we can derive from first principles. It is an emergent property of the GCMs. However it is a very robust feature of global warming simulations.
If relative humidity is nearly constant under global warming, and water vapor is a greenhouse gas, this implies a positive feedback that will amplify the warming for a given radiative forcing.
Thus far our simple models have ignored this process, and we have not been able to use them to assess the climate sensitivity.
To proceed towards more realistic models, we have two options:
We will now explore this second option, so that we can continue to think of the global energy budget under climate change as a process occurring in a single column.
We are going to adopt a parameterization first used in a very famous paper:
Manabe, S. and Wetherald, R. T. (1967). Thermal equilibrium of the atmosphere with a given distribution of relative humidity. J. Atmos. Sci., 24(3):241–259.
This paper was the first to give a really credible calculation of climate sensitivity to a doubling of CO2 by accounting for the known spectral properties of CO2 and H2O absorption, as well as the water vapor feedback!
The parameterization is very simple:
We assume that the relative humidity $r$ is a linear function of pressure $p$:
$$ r = r_s \left( \frac{p/p_s - 0.02}{1 - 0.02} \right) $$where $p_s = 1000$ hPa is the surface pressure, and $r_s$ is a prescribed surface value of relative humidity. Manabe and Wetherald set $r_s = 0.77$, but we should consider this a tunable parameter in our parameterization.
Since this formula gives a negative number above 20 hPa, we also assume that the specific humidity has a minimum value of $0.005$ g/kg (a typical stratospheric value).
This formula is implemented in climlab.radiation.ManabeWaterVapor()
Using this parameterization, the surface and tropospheric specific humidity will always increase as the temperature increases.
Here is a brief introduction to the climlab.BandRCModel process.
This is a model that divides the spectrum into 7 distinct bands: three shortwave and four longwave.
As we will see, the process works much like the familiar climlab.RadiativeConvectiveModel.
The shortwave is divided into three channels:
The longwave is divided into four bands:
The longwave decomposition is not as easily related to specific wavelengths, as in reality there is a lot of overlap between H$_2$O and CO$_2$ absorption features (as well as absorption by other greenhouse gases such as CH$_4$ and N$_2$O that we are not representing).
In [7]:
import climlab
from climlab import constants as const
First try a model with all default parameters. Usage is very similar to the familiar RadiativeConvectiveModel.
In [8]:
col1 = climlab.BandRCModel()
print( col1)
Check out the list of subprocesses.
We now have a process called H2O, in addition to things we've seen before.
The state variables are still just temperatures:
In [9]:
col1.state
Out[9]:
But the model has a new input field for specific humidity:
In [10]:
col1.q
Out[10]:
The H2O process sets the specific humidity field at every timestep to a specified profile, determined by air temperatures. More on that below. For now, let's compute a radiative equilibrium state.
In [11]:
col1.integrate_years(2)
In [12]:
# Check for energy balance
col1.ASR - col1.OLR
Out[12]:
In [13]:
fig, ax = plt.subplots()
ax.plot( col1.Tatm, col1.lev, 'c-', label='default' )
ax.plot( col1.Ts, climlab.constants.ps, 'co', markersize=16 )
ax.invert_yaxis()
ax.set_xlabel('Temperature (K)', fontsize=16)
ax.set_ylabel('Pressure (hPa)', fontsize=16 )
ax.set_title('Temperature profiles', fontsize = 18)
ax.grid()
fig
Out[13]:
By default this model has convective adjustment. We can set the adjusted lapse rate by passing a parameter when we create the model.
The model currently has no ozone (so there is no stratosphere). Not very realistic!
The Band model is aware of three different absorbing gases: O3 (ozone), CO2, and H2O (water vapor). The abundances of these gases are stored in a dictionary of arrays as follows:
In [14]:
col1.absorber_vmr
Out[14]:
Ozone and CO2 are both specified in the model. The default, as you see above, is zero ozone, and constant (well-mixed) CO2 at a volume mixing ratio of 3.8E-4 or 380 ppm.
Water vapor is handled differently: it is determined by the model at each timestep. We make the following assumptions, following a classic paper on radiative-convective equilibrium by Manabe and Wetherald (J. Atmos. Sci. 1967):
col1.relative_humidityWe need to provide some ozone data to the model in order to simulate a stratosphere. We will read in some ozone data just as we did in Lecture 8.
In [15]:
ozone = xr.open_dataset( datapath + 'som_input/ozone_1.9x2.5_L26_2000clim_c091112.nc')
In [16]:
# Take global, annual average and convert to Kelvin
weight_ozone = cos(deg2rad(ozone.lat)) / cos(deg2rad(ozone.lat)).mean(dim='lat')
O3_global = (ozone.O3 * weight_ozone).mean(dim=('lat','lon','time'))
print(O3_global)
In [17]:
fig, ax = plt.subplots()
ax.plot( O3_global*1E6, ozone.lev)
ax.invert_yaxis()
ax.set_xlabel('Ozone (ppm)', fontsize=16)
ax.set_ylabel('Pressure (hPa)', fontsize=16 )
ax.set_title('Global, annual mean ozone concentration', fontsize = 16);
fig
Out[17]:
We are going to create another instance of the model, this time using the same vertical coordinates as the ozone data.
In [18]:
# Create the column with appropriate vertical coordinate, surface albedo and convective adjustment
col2 = climlab.BandRCModel(lev=ozone.lev)
print( col2)
In [19]:
# Set the ozone mixing ratio
col2.absorber_vmr['O3'] = O3_global.values
In [20]:
# Run the model out to equilibrium!
col2.integrate_years(2.)
In [21]:
fig, ax = plt.subplots()
ax.plot( col1.Tatm, np.log(col1.lev/1000), 'c-', label='RCE' )
ax.plot( col1.Ts, 0, 'co', markersize=16 )
ax.plot(col2.Tatm, np.log(col2.lev/1000), 'r-', label='RCE O3' )
ax.plot(col2.Ts, 0, 'ro', markersize=16 )
ax.invert_yaxis()
ax.set_xlabel('Temperature (K)', fontsize=16)
ax.set_ylabel('log(Pressure)', fontsize=16 )
ax.set_title('Temperature profiles', fontsize = 18)
ax.grid(); ax.legend()
fig
Out[21]:
Once we include ozone we get a well-defined stratosphere.
Things to consider / try:
In [22]:
col3 = climlab.process_like(col2)
print( col3)
In [23]:
# Let's double CO2.
col3.absorber_vmr['CO2'] *= 2.
In [24]:
col3.compute_diagnostics()
print( 'The radiative forcing for doubling CO2 is %f W/m2.' % (col2.diagnostics['OLR'] - col3.diagnostics['OLR']))
In [25]:
col3.integrate_years(3)
In [26]:
col3.ASR - col3.OLR
Out[26]:
In [27]:
print( 'The Equilibrium Climate Sensitivity is %f K.' % (col3.Ts - col2.Ts))
In [28]:
# An example with no ozone
col4 = climlab.process_like(col1)
print( col4)
In [29]:
col4.absorber_vmr['CO2'] *= 2.
col4.compute_diagnostics()
print( 'The radiative forcing for doubling CO2 is %f W/m2.' % (col1.OLR - col4.OLR))
In [30]:
col4.integrate_years(3.)
col4.ASR - col4.OLR
Out[30]:
In [31]:
print( 'The Equilibrium Climate Sensitivity is %f K.' % (col4.Ts - col1.Ts))
Interesting that the model is MORE sensitive when ozone is set to zero.
In [32]:
%load_ext version_information
%version_information numpy, matplotlib, xarray, climlab
Out[32]:
The author of this notebook is Brian E. J. Rose, University at Albany.
It was developed in support of ATM 623: Climate Modeling, a graduate-level course in the Department of Atmospheric and Envionmental Sciences
Development of these notes and the climlab software is partially supported by the National Science Foundation under award AGS-1455071 to Brian Rose. Any opinions, findings, conclusions or recommendations expressed here are mine and do not necessarily reflect the views of the National Science Foundation.
In [ ]: